	subroutine xsplxy(xk,yk,sk,t,xt,yt,nk,nt,
     &			nna0,nfix,npa0,idadossensib)

	implicit double precision(a-h,o-z)
	INCLUDE 'param_dim.inc'
	dimension xk(nk+2),yk(nk+2),sk(nk+2),t(nt),xt(nt),yt(nt)
	integer PCnnp(nparg_max) ! conta qtos nohs parametros cada ponto de controle influencia
	integer PCnp(nparg_max,NX) ! guarda quais nohs parametros cada ponto de controle influencia
	double precision PCgama(nparg_max,NX) ! guarda os pesos com que cada coordenada dos PCs influencia cada noh parametro
	double precision PCgamaSx(nparg_max,NX),PCgamaSy(nparg_max,NX)! guarda os pesos com que cada parametro S dos PCs influencia cada Coordenada X ou Y dos nohs parametros

	common/dadossensib/PCnnp,PCnp,PCgama,PCgamaSx,PCgamaSy

	if(idadossensib .eq.1)then
		npc = nk - 1	
		do i=npa0/3+1  , npa0/3 + npc
			PCnnp(i)=0 ! zera os contadores de nohs para os PC's desta spline
		enddo

! 		do i=1,nparg_max
! 			do j=1, NX
! 				PCgamaSx(i,j) = 0.
! 				PCgamaSy(i,j) = 0.
! 			enddo
! 		enddo

	endif

	do it=1,nt	

		itk1= 2 + t(it)
		itk2= itk1+1
		
		xk0= xk(itk1-1)  ! cuidado
		xk1= xk(itk1)
		xk2= xk(itk2)

		if((itk2+1) .gt. (nk+2))then
		  xk3=xk(1)
		else
		  xk3= xk(itk2+1)  ! cuidado
		endif
		
		yk0= yk(itk1-1)  ! cuidado
		yk1= yk(itk1)
		yk2= yk(itk2)


		if((itk2+1) .gt. (nk+2))then
		  yk3=yk(1)
		else
		  yk3= yk(itk2+1)  ! cuidado
		endif		
		
		
		sk0= sk(itk1-1)  ! cuidado
		sk1= sk(itk1)
		sk2= sk(itk2)

		if((itk2+1) .gt. (nk+2))then
		  sk3=sk(1)
		else
		  sk3= sk(itk2+1)  ! cuidado
		endif	

		Tkp = itk1 + sk1
		Tk1p= itk2 + sk2
		Tk2m= itk1 - sk1
		Tk3m= itk2 - sk2
		
		pkm1= 2.D0*( (itk1-1) - Tkp )*( (itk1-1) - Tkp )
		pk  = 2.D0*(  itk1    - Tk1p)*(  itk1    - Tk1p)
		pkp1= 2.D0*(  itk2    - Tk2m)*(  itk2    - Tk2m)
		pkp2= 2.D0*( (itk2+1) - Tk3m)*( (itk2+1) - Tk3m)

		if(idadossensib .eq.1)then
			Dpkm1DTkp  = -4.0*((itk1-1) - Tkp )
			DpkDTk1p   = -4.0*( itk1    - Tk1p)
			Dpkp1DTk2m = -4.0*( itk2    - Tk2m)
			Dpkp2DTk3m = -4.0*((itk2+1) - Tk3m)
		endif

		t0= (t(it) + 2.d0 - Tkp )/( itk1-1 - Tkp )
		t1= (t(it) + 2.d0 - Tk1p)/( itk1   - Tk1p)
		t2= (t(it) + 2.d0 - Tk2m)/( itk2   - Tk2m)
		t3= (t(it) + 2.d0 - Tk3m)/( itk2+1 - Tk3m)

		if(idadossensib .eq.1)then
			Dt0DTkp  = -1.0 / (  itk1-1 - Tkp  ) + 
     &			(t(it) + 2.d0 - Tkp )/(  itk1-1 - Tkp  )**2
			Dt1DTk1p = -1.0 / (  itk1    - Tk1p) + 
     &			(t(it) + 2.d0 - Tk1p)/(  itk1    - Tk1p)**2
			Dt2DTk2m = -1.0 / (  itk2    - Tk2m) + 
     &			(t(it) + 2.d0 - Tk2m)/(  itk2    - Tk2m)**2
			Dt3DTk3m = -1.0 / ( itk2+1 - Tk3m) + 
     &			(t(it) + 2.d0 - Tk3m)/( itk2+1 - Tk3m)**2
		endif


		if(t(it)+2.d0 .gt. Tkp) then
			A0=0.d0
		else
			A0= fpt(pkm1, t0)
		endif
		A1=  fpt(pk  , t1)
		A2=  fpt(pkp1, t2)

		if(t(it)+2.d0 .lt. Tk3m) then
			A3=0.d0
		else
			A3= fpt(pkp2, t3)
		endif

		if(idadossensib .eq.1)then
			if(t(it)+2.d0 .gt. Tkp) then
				DA0Dsk1=0.d0
			else
				DA0Dsk1= +1.0*(DfptDp(pkm1, t0) * Dpkm1DTkp + 
     &					DfptDt(pkm1, t0) * Dt0DTkp)
			endif

			DA1Dsk2= +1.0*(DfptDp(pk  , t1) * DpkDTk1p + 
     &					DfptDt(pk  , t1) * Dt1DTk1p ) 
			DA2Dsk1= -1.0*(DfptDp(pkp1, t2) * Dpkp1DTk2m + 
     &					DfptDt(pkp1, t2) * Dt2DTk2m)
	
			if(t(it)+2.d0 .lt. Tk3m) then
				DA3Dsk2 =0.d0
			else
				DA3Dsk2 = -1.0*(DfptDp(pkp2, t3) * Dpkp2DTk3m + 
     &					DfptDt(pkp2, t3) * Dt3DTk3m)
			endif
		endif

		xt(it)=(A0*xk0 + A1*xk1 + A2*xk2+ A3*xk3)/(A0+A1+A2+A3)
		yt(it)=(A0*yk0 + A1*yk1 + A2*yk2+ A3*yk3)/(A0+A1+A2+A3)


		if(idadossensib .eq.1 .and. it .lt. nt)then
			! salva dados uteis para o calculo da sensibilidades
			! ip0,ip1,ip2,ip3 sao os 4 pontos de controle envolvidos
			ip3= itk1 + 2 -1 ! o -1 é para relacionar xk com a ordem dos pc dados
			if(ip3 .gt. npc) ip3 = ip3 - npc
			ip2 =ip3 - 1
			if(ip2 .lt. 1) ip2 = ip2 + npc
			ip1 =ip2 - 1
			if(ip1 .lt. 1) ip1 = ip1 + npc	
			ip0 =ip1 - 1
			if(ip0 .lt. 1) ip0 = ip0 + npc	
	! 	 	write(*,*)'it: ',it,'pontos', ip0,ip1,ip2,ip3
	
			ip0 = ip0 + npa0/3 ! passa ip local para ip global
			ip1 = ip1 + npa0/3
			ip2 = ip2 + npa0/3
			ip3 = ip3 + npa0/3
			



		
			PCnnp(ip0)=PCnnp(ip0)+1
			PCnnp(ip1)=PCnnp(ip1)+1
			PCnnp(ip2)=PCnnp(ip2)+1
			PCnnp(ip3)=PCnnp(ip3)+1

			if(npc .eq. 3)then ! neste caso, o mesmo PC foi incrementado duas vezes (no 1o e no 4o xk)
				PCnnp(ip3)=PCnnp(ip3)-1 !desfaz o incremento duplicado
			endif
		
		
			idnoh=nna0+it
! 			if(it .eq. nt)then
! 			    idnoh = nna0+1
! 			endif
			PCnp(ip0,PCnnp(ip0)) = idnoh
			PCnp(ip1,PCnnp(ip1)) = idnoh
			PCnp(ip2,PCnnp(ip2)) = idnoh
			PCnp(ip3,PCnnp(ip3)) = idnoh
		
			PCgama(ip0,PCnnp(ip0)) = A0 / (A0+A1+A2+A3)
			PCgama(ip1,PCnnp(ip1)) = A1 / (A0+A1+A2+A3)
			PCgama(ip2,PCnnp(ip2)) = A2 / (A0+A1+A2+A3)
			PCgama(ip3,PCnnp(ip3)) = A3 / (A0+A1+A2+A3)

			PCgamaSx(ip0,PCnnp(ip0)) = 0.0
			PCgamaSx(ip1,PCnnp(ip1)) = (DA0Dsk1*xk0 + DA2Dsk1*xk2)/
     &			(A0+A1+A2+A3)-(A0*xk0+A1*xk1+A2*xk2+A3*xk3)*
     &			(DA0Dsk1 + DA2Dsk1)/(A0+A1+A2+A3)**2 
			PCgamaSx(ip2,PCnnp(ip2)) = (DA1Dsk2*xk1 + DA3Dsk2*xk3)/
     &			(A0+A1+A2+A3)-(A0*xk0+A1*xk1+A2*xk2+A3*xk3)*
     &			(DA1Dsk2 + DA3Dsk2)/(A0+A1+A2+A3)**2
			PCgamaSx(ip3,PCnnp(ip3)) = 0.0	

			PCgamaSy(ip0,PCnnp(ip0)) = 0.0
			PCgamaSy(ip1,PCnnp(ip1)) = (DA0Dsk1*yk0 + DA2Dsk1*yk2)/
     &			(A0+A1+A2+A3)-(A0*yk0+A1*yk1+A2*yk2+A3*yk3)*
     &			(DA0Dsk1 + DA2Dsk1)/(A0+A1+A2+A3)**2 
			PCgamaSy(ip2,PCnnp(ip2)) = (DA1Dsk2*yk1 + DA3Dsk2*yk3)/
     &			(A0+A1+A2+A3)-(A0*yk0+A1*yk1+A2*yk2+A3*yk3)*
     &			(DA1Dsk2 + DA3Dsk2)/(A0+A1+A2+A3)**2
			PCgamaSy(ip3,PCnnp(ip3)) = 0.0	

! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! 
! ! ! ! !  Verificacao por diferencas finitas
! ! ! ! 			deltaS=0.01
! ! ! ! 	
! ! ! ! 			do i=0,3
! ! ! ! 
! ! ! ! 				if(i.eq.0)then
! ! ! ! 					sk0 = sk0 + deltaS
! ! ! ! 					ip = ip0
! ! ! ! 				else if(i.eq.1)then
! ! ! ! 					sk0 = sk0 - deltaS
! ! ! ! 					sk1 = sk1 + deltaS
! ! ! ! 					ip = ip1
! ! ! ! 				else if(i.eq.2)then
! ! ! ! 					sk1 = sk1 - deltaS
! ! ! ! 					sk2 = sk2 + deltaS
! ! ! ! 					ip = ip2
! ! ! ! 				else if(i.eq.3)then
! ! ! ! 					sk2 = sk2 - deltaS
! ! ! ! 					sk3 = sk3 + deltaS
! ! ! ! 					ip = ip3
! ! ! ! 				endif
! ! ! ! 
! ! ! ! 				Tkp = itk1 + sk1
! ! ! ! 				Tk1p= itk2 + sk2
! ! ! ! 				Tk2m= itk1 - sk1
! ! ! ! 				Tk3m= itk2 - sk2
! ! ! ! 			
! ! ! ! 				pkm1= 2.D0*( (itk1-1) - Tkp )*( (itk1-1) - Tkp )
! ! ! ! 				pk  = 2.D0*(  itk1    - Tk1p)*(  itk1    - Tk1p)
! ! ! ! 				pkp1= 2.D0*(  itk2    - Tk2m)*(  itk2    - Tk2m)
! ! ! ! 				pkp2= 2.D0*( (itk2+1) - Tk3m)*( (itk2+1) - Tk3m)
! ! ! ! 	
! ! ! ! 				t0= (t(it) + 2.d0 - Tkp )/( itk1-1 - Tkp )
! ! ! ! 				t1= (t(it) + 2.d0 - Tk1p)/( itk1   - Tk1p)
! ! ! ! 				t2= (t(it) + 2.d0 - Tk2m)/( itk2   - Tk2m)
! ! ! ! 				t3= (t(it) + 2.d0 - Tk3m)/( itk2+1 - Tk3m)
! ! ! ! 
! ! ! ! 				if(t(it)+2.d0 .gt. Tkp) then
! ! ! ! 					A0=0.d0
! ! ! ! 				else
! ! ! ! 					A0= fpt(pkm1, t0)
! ! ! ! 				endif
! ! ! ! 					A1=  fpt(pk  , t1)
! ! ! ! 					A2=  fpt(pkp1, t2)
! ! ! ! 		
! ! ! ! 				if(t(it)+2.d0 .lt. Tk3m) then
! ! ! ! 					A3=0.d0
! ! ! ! 				else
! ! ! ! 					A3= fpt(pkp2, t3)
! ! ! ! 				endif
! ! ! ! 
! ! ! ! 				xtaux=(A0*xk0 + A1*xk1 + A2*xk2+ A3*xk3)/(A0+A1+A2+A3)
! ! ! ! 				ytaux=(A0*yk0 + A1*yk1 + A2*yk2+ A3*yk3)/(A0+A1+A2+A3)
! ! ! ! 
! ! ! ! 				PCgamaSx(ip,PCnnp(ip)) = PCgamaSx(ip,PCnnp(ip))+
! ! ! !      &					(xtaux - xt(it))/ deltaS
! ! ! ! 				PCgamaSy(ip,PCnnp(ip)) = PCgamaSy(ip,PCnnp(ip))+
! ! ! !      &					(ytaux - yt(it))/ deltaS
! ! ! ! 			enddo
! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! 
		endif
	
	end do

! 	if(idadossensib .eq.1)then
! 		write(*,*)'xsplxy - PCGama'
! 		do i=1,npc
! 			write(*,*)(PCGama(i,j),j=1,PCnnp(i))
! 		enddo
! 	endif
	
	return
	end













      double precision function fpt(p,t)
      implicit double precision(a-h,o-z)
         
          fpt = (10.d0 - p + (2.d0*p-15.d0)*t + (6.d0-p)*t*t )*t*t*t
       return
       end

      double precision function DfptDp(p,t)
      implicit double precision(a-h,o-z)
         
          DfptDp = -t**3 + 2.0*t**4 - t**5

       return
       end

      double precision function DfptDt(p,t)
      implicit double precision(a-h,o-z)
         
          DfptDt = 30.0*t**2 - 3.0*p*t**2 +
     *             4.*(2.*p -15.)*t**3 + 5.0*(6.0 - p)*t**4

       return
       end